Phonon spectrum, thermodynamic properties, and pressure-temperature phase 

diagram of uranium dioxide 
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We present a study of the structural phase transition, mechanical and thermodynamic properties 
of UO2 by means of the local density approximation (LDA)+?7 approach. A phase transition 
pressure of 40 GPa, which agrees well with the experimental value of 42 GPa, is obtained from 
theory at K. Pressure-induced enhancements of elastic constants, elastic moduli, elastic wave 
velosities, and Debye temperature of the fluorite phase are observed. Phonon spectrums of both 
the ground state fluorite structure and high pressure cotunnite structure calculated by the supercell 
approach show that the cotunnite structure is dynamically unstable under ambient pressure. Based 
on the imaginary mode along the F—X direction and soft phonon mode along the T—Z direction, a 
transition path from cotunnite to fluorite has been identified. We calculate the lattice vibrational 
energy in the quasiharmonic approximation using both first-principles phonon density of state and 
the Debye model. Calculated temperature dependence of lattice parameter, entropy, and specific 
heat agree well with experimental observations in the low temperature domain. The difference 
of Gibbs free energy between the two phases of UO2 has predicted a boundary in the pressure- 
temperature phase diagram. The solid-liquid boundary is approximated by an empirical equation 
using our calculated elastic constants. 

PACS numbers: 71.27.+a, 61.50.Ks, 62.20.-x, 63.20.dk 



I. INTRODUCTION 



Due to its critical importance in nuclear fuel cycle and 
complex electronic structure arisen from a partially oc- 
cupied 5/ orbital, uranium dioxide (UO2) has been stud- 
ied extensively in experiments [__-_5j and computational 
simulations [6h14|. The 5/ electrons in UO2 play a piv- 
otal role in understanding its electronic, thermodynamic, 
and magnetic properties [l5j |. Using conventional den- 
sity functional theory (DFT) of exchange-correlation po- 
tential, i.e., the local density approximation (LDA) or 
generalized gradient approximation (GGA), an incorrect 
ferromagnetic (FM) conducting ground state of UO2 was 
observed Q since an error produced by underestimat- 
ing the strong on-site Coulomb repulsion of the 5/ elec- 
trons. Same problems have been confirmed in previous 
investigations of NpC>2 [16] and PuC>2 [H[ within the 
pure LDA/GGA schemes. Fortunately, for Pu02 a the- 
ory based on completely localized 5/ states reproduced 
well the crystal field splittings as well as the magnetic 
susceptibility [17(. The f—>f antiferromagnetic (AFM) 
Mott-Hubbard insulator nature of UO2 has been well 
reproduced in LDA/GGA+t/ Q, hybrid density func- 
tional of (Heyd, Scuseria, and Enzerhof) HSE self- 
interaction corrected local spin-density (SIC-LSD) [Hj], 
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and LDA plus Dynamical Mean-Field Theory (DMFT) 
[ill calculations to account for the photoelectron spec- 
troscopy experiments [H, __| . 

At ambient conditions, UO2 crystallizes in a cubic flu- 
orite structure (Fm3m, No. 225) with cations arranging 
in a face-centered cubic (fee) structure and anions occu- 
pying tetrahedral sites. Similar to the high-pressure be- 
havior of TI1O2 and Pu02 [l9j, a recent hydrostatic com- 
pression experiment [3J has shown that UO2 also trans- 
forms to the orthorhombic structure of cotunnite-type 
(Prima, No. 62) at room temperature, beyond 42 GPa. 
This kind of pressure-induced phase transition (PT) for 
actinidc dioxides is the same as for the alkaline earth 
fluorides [20j and has not been sufficiently studied, al- 
though many experiments [!, [l9[ and theoretical works 
[I3L I2H l22l | have paid attentions on this issue. The co- 
tunnite phase data are scarce in the literature, especially 
for its thermodynamic properties and vibrational char- 
acters. Temperature contributions to the PT have not 
been included in previous studies. On the other hand, 
the melting performances of UO2 also have not been well 
investigated. Only few experiments have been conducted 
to describe the melting features of UO2 near ambient 
pressure, because of the difficult experimental conditions 
required to control and monitor the PT [231 ] . 

In our previous systematic work jl3| . we have pre- 
sented structural, electronic, and mechanical properties 
of AFM UO2 in its ground-state fluorite phase and high- 
pressure cotunnite phase at their corresponding equilib- 
rium states within LDA+ U formalism with U —A eV. 
The lattice parameter ao=5.449 A and bulk modulus 



2 



B=220.0 GPa for Fm3m U0 2 obtained by the third- 
order Birch-Murnaghan equation of state (EOS) [HJ fit- 
ting were found to be in perfect agreement with results of 
recent LDA+ U calculation [2I (a =5.448 A and 5=218 
GPa) and experiments 0, [26| (ao=5.47 A and 5=207 
GPa). In the present work, we perform an extended 
study of the structural, mechanical, and thermodynamic 
properties of UO2 in the pressure range from to 250 
GPa and in a temperature interval from to 4000 K by 
employing the LDA+ U and GGA+ U schemes as imple- 
mented by Dudarev et al. @, US HH . The total energies 
of nonmagnetic (NM), AFM, and FM phases of the fluo- 
rite structure have been calculated in a wide range of the 
effective Hubbard U parameter within LDA/GGA+[7 
schemes to guarantee the validity of the ground-state cal- 
culations. At K, a Fm3m—i'Pnma PT pressure of 40 
GPa is predicted. So, we have calculated the elastic con- 
stants, elastic moduli, Poisson's ratio, elastic wave ve- 
locities, and Debye temperature of AFM fluorite UO2 in 
the pressure range from to 40 GPa. The structural 
transition path of cotunnite phase to fluorite phase as 
well as the melting behavior have been studied based 
upon our calculated phonon dispersions, Gibbs free en- 
ergy, and elastic constants. Thermodynamic properties 
including of Gibbs free energy, temperature dependence 
of the lattice parameter and the bulk modulus, entropy, 
and specific heat have also been evaluated. The rest of 
this paper is arranged as follows. In Sec. II the compu- 
tational methods are described. In Sec. Ill we present 
and discuss our results. In Sec. IV we summarize the 
conclusions of this work. 



II. COMPUTATIONAL METHODS 
A. Computational details 

First-principles DFT calculations on the basis of the 
frozen-core projected augmented wave (PAW) method of 
Blochl [2^ are performed within the Vienna ab initio 
simulation package (VASP) [30| . where the exchange and 
correlation effects are described by the LDA and GGA 
[3l|, H2]. For the plane- wave set, a cutoff energy of 500 eV 
is used. The fc-point meshes in the full wedge of the Bril- 
louin zone (BZ) are sampled by 9 x 9 x 9 and 9x15x9 grids 
according to the Monkhorst-Pack (MP) [33| scheme for 
fluorite and cotunnite UO2, respectively, and all atoms 
are fully relaxed until the Hellmann-Feynman (HF) forces 
become less than 0.02 eV/A. The U 6s 2 7s 2 6p 6 6d 2 5f 2 and 
the O 2s 2 2p 4 orbitals are treated as valence electrons. 
Similar to our previous studies [HEl], 

the strong on-site 

Coulomb repulsion among the localized U 5/ electrons is 
described by using the LDA/GGA+ U formalisms for- 
mulated by Dudarev et al. 0, H3, HH , where the double 
counting correction has already been included. In this 
paper the Coulomb U is treated as a variable, while the 
exchange energy is set to be a constant J=0.51 eV. This 
value of J had been checked carefully by Dudarev et al. 
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FIG. 1: (Color online) Dependence of the total energies (per 
formula unit) on U for NM, FM, and AFM UO2 in Fm3m 
phase. 



0, H3, [H- Since only the difference between U and J 
is significant (28[, we will henceforth label them as one 
single parameter, for simplicity labeled as U, while keep- 
ing in mind that the non-zero J has been used during 
calculations. 



Both spin-unpolarized and spin-polarized calculations 
are performed in this study. We show in Fig. 1 the de- 
pendence of the total energy (per formula unit at respec- 
tive optimum geometries) on U for NM, FM, and AFM 
phases within the LDA+ U and GGA+ U formalisms. 
Clearly, the NM phase is not energetically favorable in 
the LDA+ U and GGA+ U formalisms, compared to FM 
and AFM phases. Therefore, the results of NM calcu- 
lations are not presented in the following. At [7=0 and 
1.0 eV, the total energy of the FM phase is lower than or 
almost equal to that of the AFM phase either in LDA+ U 
scheme or GGA+ U scheme. However, as shown in Fig. 
1 , it is clear that the total energy of the AFM phase de- 
creases to become lower than that of the FM phase when 
increasing U. The total-energy differences (Efm—Eafm) 
within the LDA+U and GGA+U at U=A eV are 0.755 
and 1.053 eV, respectively. The AFM stable nature es- 
tablished by LDA+ U or GGA+ U scheme is consistent 
with experiments 0, H3, HH . In the following study, we 
focus on results of the AFM phase. Also, although the 
spin-orbit coupling (SOC) is important for certain prop- 
erties of heavy metal compounds, inclusion of SOC has 
limited effect on the bulk-properties and chemical bind- 
ing of UO2 and Pu02 [8j, [l3j, I3q - l39j . The reason for this 
is that the /-shell is localized in these compounds, and 
the chemical binding is provided by the spd-st&tes of U 
and the sp-states of O, and for these states SOC is less 
important. Therefore, in our present work of U0 2 , the 
SOC is not included. 



3 



B. Elastic properties, Debye temperature, and 
melting temperature 

To avoid the Pulay stress problem, the geometry op- 
timization at each volume is performed at fixed volume 
rather than constant pressure. Elastic constants for cubic 
symmetry (Cn, C12, and C44) and orthorhombic struc- 
ture (Cn, C12, C13, C 22 , G23, C 33 , C44, C 55 , and C 66 ) are 
calculated by applying stress tensors with various small 
strains onto the equilibrium structures. The strain ampli- 
tude S is varied in steps of 0.006 from S=— 0.036 to 0.036. 
Detailed descriptions of calculation scheme please see our 
previous work After obtaining elastic constants, the 
polycrystallinc bulk modulus B and shear modulus G 
are calculated from the Voigt-Reuss-Hill (VRH) approx- 
imations [iof . The Young's modulus E and Poisson's 
ratio v are calculated through E — 9BG/(3B + G) and 
v = (3B - 2G)/[2(3B + G)]. In calculation of the Debye 
temperature (0d), we use the relation 
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k B J Vm ' 



(1) 



where h and k B are Planck and Boltzmann constants, 
respectively, n is the number of atoms in the molecule, f2 
is molecular volume, and v m is the average sound wave 
velocity. The average wave velocity in the polycrystallinc 
materials is approximately given as 
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1-1/3 



(2) 



where vt = y/G/p (p is the density) and v\ = 
y (3B + 4G)/3p are the transverse and longitudinal elas- 
tic wave velocity of the polycrystalline materials, respec- 
tively. The melting temperature (T rn ) in units of K for 
cubic U0 2 are deduced from elastic constant (Cn) by an 
approximate empirical linear formula fill ]: 



T m = 553 + 5.91Cn, 



(3) 



where the Cn is in units of GPa and the standard error 
is about ±300 K. 



C. Phonon and thermodynamic properties 

We use the supercell approach (42[ and the small dis- 
placement method as implemented in the FROPHO code 
[43j to calculate the phonon curves in the BZ and the cor- 
responding phonon density of states (DOS) for both flu- 
orite and cotunnite phases of UO2 . In the interpolation 
of the force constants for the phonon dispersion curve 
calculations, 3x3x3 and 3x5x3 MP fc-point meshes are 
used for Fm3m 2x2x2 and Prima 2x2x2 supercells, re- 
spectively. The forces induced by small displacements 
are calculated within VASP. 



Thermodynamic properties can be determined by 
phonon calculation using the quasiharmonic approxima- 
tion (QHA) [13, 13] or by the quasiharmonic Debye model 
45]. Within these two models, the Gibbs free energy 
G(T,P) is written as 



G(T, P) = F(T, V) + PV. 



(4) 



Here, F(T,V) is the Helmholtz free energy at tempera- 
ture T and volume V and can be expressed as 

F(T,V) = E(V) + F mb (T, V) + F el (T, V), (5) 

where E(V) is the ground-state total energy, F V i b (T, V) is 
the vibrational energy of the lattice ions and F e i ( T, V) is 
the thermal electronic contribution. Similar to our pre- 
vious work 13], F e i(T, V) is also not included in present 
study. 

Under QHA, the F v n,(T,V) can be calculated by 



F vib {T,V) = k B T f 
Jo 



g(iv) In 



2 sinh 



2k B T 



(6) 

where to represents the phonon frequencies and g{uj) is 
the phonon DOS. This formula requests positive results 
of the phonon DOS. So it is not suitable for dynamically 
unstable phases. Instead, the vibration energy for phases 
with imaginary phonon frequencies can be estimated by 
the Debye model 



F mb {T,V) = h B 6 D +k B T 



3 In 1 - 



T 



(7) 

where ^k B 8 B is zero-point energy due to lattice ion vi- 
bration at K and D(6 B) /T) the Debye integral written 

as D{9 D /T) = 3/(9 D /T) 3 f° D/T x 3 /(e x - l)dx. For a 
more detailed computation scheme of the Debye model 
please see Ref. [45j . 



III. RESULTS 

A. Phase transition at K 

In the present study, the structural parameter and bulk 
modulus for Frr&m AFM U0 2 calculated using LDA+ U 
formalism with U =4 eV are identical to our previous 
results [131 ]. The insulating energy band gap (E g ) and 
spin moment per U atom {p, mag .) are calculated to be 1.9 
eV and 1.977 p B , respectively, which are in g ood agree- 
ment with previous LDA+U calculation [21] (E g =1.45 
eV and /i ma g. = l-93 ur) an d experiments (E g =2.0 eV 
and ^ mQ5 .=1.74 p B [35]). For Pnma U0 2 in AFM phase, 
we obtain the optimized structural lattice parameters a, 
6, and c to be 5.974, 3.604, and 6.967 A, respectively, 
after fitting the energy-volume data to the EOS. Band 
gap and spin moment are calculated to be 1.6 eV and 
2.002 p B , respectively. Thus, the band gaps should not 
increase from Fm3m phase to Pnma phase when using 
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TABLE I: Lattice constants, elastic constants, various moduli, Poisson's ratio (v), density (p), transverse (ft), longitudinal (vi) 
and average (v m ) sound velocities, and Debye temperature (6d) for FmZm AFM UO2 at different pressures calculated within 
LDA+ U formalism with U =4 eV. For comparison, experimental values and other calculation results at GPa are also listed. 



Pressure 


a 


C11 


c 


C 


B 


G 


E 


V 




r 




V) 


Vm 


0d 


(GPa) 


(A) 


(GPa) 


(GPa) 


(GPa) 


(GPa) 


(GPa) 


(GPa) 




( 2/cm^) 
v&/ J 


( m/s 1 ) 


fm/s) 


fm/s) 


(K) 





5.449 


389.3 


138.9 


71.3 


222.4 


89.5 


236.8 


0.323 


11.084 


2841.8 


5552.7 


3183.4 


398.1 


5 


5.408 


414.8 


154.5 


94.3 


241.2 


107.3 


280.4 


0.306 


11.343 


3076.1 


5821.1 


3438.7 


433.3 


10 


5.373 


438.2 


166.7 


106.7 


257.2 


117.5 


305.9 


0.302 


11.565 


3187.5 


5982 3 


3561.2 


451.7 


15 


5.340 


459.2 


181.6 


118.9 


274.1 


126.5 


328.9 


0.300 


11.776 


3277.7 


6131.8 


3661.1 


467.2 


20 


5.310 


479.8 


195.6 


131.3 


290.4 


135.5 


351.8 


0.298 


11.979 


3363.3 


6270.7 


3755.8 


482.0 


25 


5.282 


500.3 


208.3 


143.6 


305.6 


144.6 


374.6 


0.296 


12.175 


3445.7 


6398.0 


3846.7 


496.3 


30 


5.254 


520.8 


221.6 


156.0 


321.3 


153.4 


397.0 


0.294 


12.364 


3522.3 


6521.7 


3931.4 


509.8 


35 


5.229 


540.0 


233.7 


167.8 


335.8 


161.8 


418.2 


0.292 


12.546 


3591.1 


6630.0 


4007.3 


522.2 


40 


5.204 


558.0 


246.7 


180.9 


350.4 


170.3 


439.8 


0.291 


12.722 


3659.1 


6737.7 


4082.5 


534.5 


Expt. 


5.4731 a 


389.3 6 


118. 7 b 


59.7* 


209.0 i> 


83.0 6 


221. 6 


0.324'' 










385 6 , 395 c 


LDA+ U d 


5.448 


380.9 


140.4 


63.2 


220.6 


82.0 


218.9 


0.335 










399 



" Reference 0], b Reference [H, c Reference [H], d Reference 



same value of Hubbard parameters. This is different from 
previous LDA+ U calculation [2l| , where an increase of 
the band gap was found at a cell volume close to the 
transition pressure from 0.8 eV in the fluorite phase to 
2.4 eV in the cotunnite phase, by using different values of 
Hubbard parameters. Although the PT behaviors should 
seek help from larger value of Hubbard U parameter (see 
follows), we present equilibrium state properties of Prima 
U0 2 with [7=4 eV. 

In Fig. 2, we show the total energy vs cell volume 
curves since their critical role in solid-state behavior de- 
scriptions. Because the pressure induced PT can not be 
properly modeled by [/=4 eV within LDA+ U, we calcu- 
late the energy- volume data of Pnma phase by [/=5.5 
eV. Detailed explanation on this issue can be found in 
Ref. [2l|. It is clear that the Fm3m phase is stable un- 
der ambient condition and upon compression it will tran- 
sit to the Pnma phase. At K, the Gibbs free energy is 
equal to the enthalpy H. After calculation, we can plot in 
Fig. 2 (see the inset) the relative enthalpies of the Pnma 
phase with respect to the Fm3m phase as a function of 
pressure. The cross point predict a PT pressure of 40 
GPa, which is consistent with previous LDA+ U calcula- 
tion [2l| value of ~38 GPa and also the experimentally 
observed 42 GPa [|. 



B. Elasticity of fluorite UO2 

Elastic constants can measure the resistance and me- 
chanical features of crystal to external stress or pres- 
sure, thus describing the stability of crystals against 
elastic deformation. We present in Table I the lattice 
constants, elastic constants, bulk moduli, shear moduli, 
Young moduli, Poisson's ratio, density, elastic wave ve- 
locities, and Debye temperature for Fm3m AFM UO2 
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FIG. 2: (Color online) Total energy vs the cell volume for 
AFM UO2 in FmZrn and Pnma phases. Results of FmZm are 
calculated within LDA+ U formalism with U =4 eV, while for 
Pnma phase are obtained by U =5.5 eV. A PT at 40 GPa is 
predicted by the pressure dependence of the enthalpy differ- 
ences of Pnma phase with respect to Fm3m phase, as shown 
in the inset. 



at different pressures. All these results are calculated 
within the LDA+ U formalism with U—4 eV. The val- 
ues at zero pressure is the same as in our previous work 
[HI]. Actually, elastic constants at GPa have been 
widely studied either by experiments [46| or by calcu- 
lations [H, H3, |48| . Our calculated results at GPa are 
consistent with corresponding values from experiments 
[4(| and recent LDA+ U work [3^] , where Sanati et al. 
have presented carefully discussions on the zero pressure 
results. In this study, various moduli, Poisson's ratio (v ), 
density (p), elastic wave velocities (ft, Vi, and v m ), and 
Debye temperature (#d) of low-pressure fluorite phase 
UO2 in the pressure range from GPa to 40 GPa have 
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been deduced from elastic constants. In the entire pres- 
sure range, Cn is prominently larger than C12, indicat- 
ing that the bonding strength along the [100]/ [010]/ [001] 
directions is clearly stronger than that of the bonding 
along the [011]/[101]/[110] directions. In fact, there are 
eight U— O covalent bonds per formula unit for fluorite 
UO2. The angle of all eight bonds with respect to the 
[100]/[010]/[001] directions is 45°. However, only four 
bonds make an angle of 45° with the [011]/ [101]/ [110] 
directions. Four other bonds are vertical to the strain 
directions of [011]/[101]/[110]. The bonds vertical to 
the strain directions have no contributions on the elastic 
strength. Therefore, the fact that C\\ > C12 for cubic 
UO2 is understandable. This kind of bonding analysis 
has been used previously to explain the different theo- 
retical tensile strengths in the three typical crystalline 
orientations of PuC>2 [Hj]. For the Debye temperature, 
our calculated result of 398.1 K is in excellent agreement 
with experiments [46|, |49( . 

As indicated in Table I, except poisson's ratio, 
pressure-induced enhancements of elastic constants, elas- 
tic moduli, elastic wave velocities, and Debye tempera- 
tures are evident. They all increase linearly with pres- 
sure. While C12 and C44 have the same increase rate of 
~2.7, Cn has a larger one of ~4.2. This also can be un- 
derstood from the previous bonding analysis. The rates 
with which B, G, and E increase, are 3.2, 2.0, and 5.1, 
respectively. Considering By = Br — (Cn + 2Ci2)/3, 
G v = (C n - C12 + 3C 44 )/5, and G R = 5(C U - 
Ci2)C44/[4C44 + 3(Cn — C12)] for cubic symmetry, we 
can understand the fact that the increase rate of G is 
only about 60% of the increase rates of B. For w t and vi , 
increase rates of 20.4 and 29.6 m/s/GPa are obtained, 
respectively. Relatively larger increase rate of transverse 
sound velocity upon compression is due to the larger en- 
hancement of the bulk modulus B with respect to the 
shear modulus G. The pressure-dependent linearly in- 
creasing behavior of Debye temperature is also clear. 
This kind of analysis for Debye temperature should sup- 
ply useful informations in practical application and/or 
theoretical analysis for UO2. 




FIG. 3: (Color online) Phonon dispersion curves (left panel) 
and corresponding phonon DOS (right panel) for UO2 in (a) 
Fm3m phase and (b) Pnma phase. All results are calculated 
within LDA+ U formalism with U =4 eV. The solid lines are 
calculated without including the polarization effects, while 
the dashed lines including. The hollow circles present the 
experimental data from Ref. [4!|. 
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C. Phonon dispersion 

The calculated phonon dispersion curves as well as 
corresponding phonon DOSs are displayed in Fig. [3] 
for Fm3m and Pnma UO2 in the AFM configura- 
tion. To our knowledge, no experimental or theoret- 
ical phonon frequency results have been published for 
the high-pressure phase of actinide dioxides. For Fm3m 
UO2, inelastic neutron scattering (49|, infrared and Ra- 
man spectroscopic [l], l50l. 51 1 experiments as well as 
LDA+DMFT [u[, MD gjand LDA/GGA+J7+SOC 



[39j calculations have been performed to picture its vibra- 
tional features. In addition to their phonon curves along 
<001>, <110>, and <111> directions in the BZ, we have 
shown phonon dispersions along T—X—K—T—L—X—W 



FIG. 4: (Color online) Schematic illustrations of the struc- 
tural transition from Pnma phase to Fm3m structure. For 
clarity, only uranium atoms are presented and atoms within 
the Pnma unit cell are labeled by star symbols. 



directions. The r~ X, T—K, and T—L lines are along 
<001>, <110>, and <111> directions, respectively. As 
having been demonstrated by Sanati et al. (39|, includ- 
ing of SOC dose not significantly modify the phonon 
DOSs and the contribution of the lattice distortion is 
also limited, however, when Hubbard corrections are not 
included, the higher-frequency phonon modes are shifted 
to lower frequencies and the optical modes are under- 
estimated. Therefore, we have calculated phonon dis- 



6 



-124 
-126 
-128 

> 

03 

— -130 
>. 

2 -132 

0) 
0) 

£ -134 
-136 



-138 



5.2 




5.3 



5.4 



5.5 
a (A) 



5.6 



5.7 



5.8 



FIG. 5: Dependence of the free energy F(T,V) on crystal 
lattice parameter a for a number of selected temperatures for 
AFM U0 2 calculated within LDA+U formalism with [7=4 
eV. 
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persions and phonon DOSs of cubic UO2 by LDA+U 
scheme without including SOC. From Fig. Eta), one can 
find that including polarization effect is necessary to cor- 
rectly account for the LO-TO splitting near the T point 
in BZ. Here, the Born effective charges (Z{j=5.54 and 
Z5=-2.77) of U and O ions for Fm3m AFM U0 2 are 
also calculated. Our phonon dispersions are overall in 
good agreement with the inelastic neutron scattering ex- 
periment [!§] and previous calculations [ill US • 



As for Pnma UO2 , we have shown phonon dispersions 
along T-X-U-Z-T-Y-T-R directions. The nota- 
tions of the high-symmetry points are T (0, 0, 0), X (0, 
i, 0), U (0, i, i), Z (0, 0, i), Y (-i, 0, 0), T (-i, 
0, i), and R (-i, |, i). Although mechanically stable 
nature of Pnma UO2 at its equilibrium state has been 
predicted by the elastic constants in our previous work 
(13| . Fig. EJb) clearly shows that the transverse acous- 
tic (TA) mode close to T point becomes imaginary along 
the T—X (i.e., the <010>) direction. This means that 
the high-pressure phase of UO2 is dynamically unstable 
under ambient pressure. In addition, we can find a clear 
soft phonon mode along T~Z (i.e., the <001>) direc- 
tion. Thus, uranium atoms in Pnma structure are easy 
to move along <010> and <001> directions. Based upon 
these observations, we show in Fig. @]a suggested path of 
the Pnma^FmSm transition. The Pnma phase can be 
viewed as an AB periodically layered structure along the 
[100] direction. In transition, firstly the adjacent (100) 
planes slip relatively along the [001] direction to create 
a face centered orthorhombic structure (as indicated by 
the arrows in Fig. @|, and second, the cell expands along 
the [010] direction and shrinks in the vertical directions 
to form the fee fluorite structure. 



FIG. 6: (Color online) Temperature dependences of (a) lattice 
parameter a(T) and (b) bulk modulus B(T) of UO2. Exper- 
imental results from [53 ] and [54|] as well as the MD results 
from [55| are also shown in panel (a). 



D. Thermodynamic properties and P — T phase 
diagram 

Calculated free energy F(T, V) curves of UO2 for tem- 
peratures ranging from up to 2000 K are shown in 
Fig. El Note that in the calculation of F(T,V), the 
ground-state total energy and phonon free energy should 
be calculated by constructing several 2x2x2 fee super- 
cells. This kind of calculation is computationally very 
expensive. In Fig. [SJ the equilibrium lattice parame- 
ters at different temperature T are also presented. The 
equilibrium volume V(T) and the bulk modulus B(T) are 
obtained by EOS fitting. Figure [5] shows the temperature 
dependence of the lattice parameter and the bulk modu- 
lus. Experimental results from [HH and [HH as well as the 
MD results from [55[ are also plotted. We observe good 
agreement of calculated lattice parameters with respect 
to the experiments in the low temperature domain, but 
somewhat lower values compared to experiments in tem- 
perature higher than 800 K. The differences may come 
from the thermal electronic contribution and/or anhar- 
monic effects. For the bulk modulus B(T), similar as 
for Pu02 [l3j , a decreasing behavior upon elevating tem- 
perature is observed. The decreasing amplitude of B(T) 
from K to 1500 K for U0 2 is -26.8 GPa, which is larger 
than that of Pu02 by about 6.2 GPa. This means that 
UO2 will be softened quicker upon increasing tempera- 
ture in comparison with Pu02- 

Using QHA and the Debye model, we have calculated 
the Gibbs free energy (G), entropy (S), and specific heat 
at constant volume (Cy) for FmZm and Pnma phases of 
UO2. Since the specific heat at constant pressure (Cp) 
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FIG. 8: (Color online) (a) Temperature dependences of Gibbs free energy curves for UO2 in Fm3m and Pnma phases at 
GPa. (b) The difference of the Gibbs free energies of the Fm3m and Pnma phases of UO2 as a function of pressure for several 
temperatures (in K). A positive value indicates that the Pnma phase is more stable. 



is similar to Cy [13|, we will not present Cp results in 
this work. Entropy and heat capacity at constant vol- 
ume results at GPa are presented in Fig. [7] and Gibbs 
free en ergy results are plotted in Fig. [H Experimental 
results [Fy, [57} of S and Cy are also showed for compar- 
ison. Under QHA, the entropy S curves for Fm3m and 
Pnma phases are almost identical to each other. Sim- 
ilar to previous calculation |39|, the S of fluorite UO2 
is underestimated in a wide range of temperatures with 
respect to experiments. However, as clearly indicated in 
Fig. Ufa), the Debye model can give proper results for 
Fm3m U02- Using the Debye model, the S curves for 
Fm3m and Pnma phases will separate along with in- 
creasing temperature. The difference between QHA and 
Debye model is due to the fact that QHA uses harmonic 
approximation and the Debye model has considered some 
anharmonic contributions in calculations of S and Cy. 
Although the Debye model is less accurate, it can supply 
a qualitative picture or even quantitative description of 
the thermodynamic properties. As shown in Fig. [7jb), 



the Cy of Fm3m UO2 under QHA agrees well with ex- 
periments up to room temperature and becomes close to 
a constant of the Dulong-Petit limit [58[ . Same trends of 
Cy curve for FmZm phase has been observed by Sanati 
et al. J39( and our results based on phonon DOS indi- 
cate that the Cy curves for Fmim and Pnma phases 
are almost identical to each other. However, under the 
Debye model, slower increasing behavior of Fm3m phase 
Cy with increasing temperature is observed with respect 
to the Pnma phase. The Debye model gives #£>=390.6 
and 352.8 K for FmZm and Pnma phases, respectively, 
which are in good agreement with 398.1 and 343.7 K [131 ] 
computed by elastic constants. 

As shown in Fig. EJa), the cross between Fm3m and 
Pnma phases Gibbs free energy curves of UO2 by Debye 
calculations, readily gives a FmSm^-Pnma phase tran- 
sition temperature of 2069 K. This implies a significant 
temperature contribution for the FmZm^Pnma phase 
transition, which is not only pressure driven. To predict 
the phase boundary of FmSm^-Pnma, we calculate here 
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the Gibbs free energy of the Fm3m and Pnma crystal 
structures in a temperature range from to 3000 K, and 
the pressure effect is studied in the range of to 45 GPa. 
The Gibbs energy difference between the fluorite and co- 
tunnitc structures (AG) of UO2 as a function of pressure 
for several temperatures is reported in Fig. Efb). At 
K, the Fm3m—>Pnma PT pressure is predicted to be 
40 GPa, which is identical to our aforementioned result. 
Along with increasing temperature in the range from 
to 2069 K, the pressure of the Fmim^Pnma transition 
decreases slightly. At higher temperatures, the Fm3m 
phase is only stable in middle pressure range. 

Once the free energies of the two experimentally ob- 
served structures are determined, the phase boundary 
can be obtained by equating Gibbs free energies at a 
given pressure and temperature. Besides, the solid-liquid 
boundary can be featured by the melting temperature 
T m . Based on these results, we can plot in Fig. [5] the 
phase diagram of UO2. Other theoretical pll. l55j and 
experimental [!, [23| values are also presented for com- 
parison. For the phase boundary between Fm3m and 
Pnma, only one point at ambient condition was reported 
in experiment. Our predicted results call for further ex- 
perimental and theoretical works, to investigate the ac- 
curacy of the theory. For the solid-liquid boundary, the 
experimentally determined melting value at zero pressure 
is 3147 K [23j. We note that experiment and recent MD 
calculation have reported the relationship between melt- 
ing point and pressure to be T m .p =3147+92. 9P and 
T m ,p =3178+115P, respectively, where P is pressure in 
unit of GPa. Using our calculated data in Fig. |H1 we 
obtain T m ^ P =2882+24. 8P. The melting point at zero 
pressure is underestimated by about 265 K, which is the 
same with previous LDA+C/ calculation [39] . The in- 
creasing rate of T m ^p is largely underestimated compared 
to experiment. Although we do not wish to state that our 
calculations are more reliable than experiment, it stands 
clear that since the experiment was performed in a nar- 
row range of pressure between 0.01 and 0.25 GPa, the 
differences between our calculation and previous experi- 
ment call for more theoretical and experimental works. 



model. Pressure-dependent linearly increasing behaviors 
of mechanical properties and Debye temperature of fluo- 

4000 n 1 1 1 1 1 1 1 , — 




FIG. 9: (Color online) Calculated P-T phase diagram for 
uranium dioxide compared with selected experimental points 
(square for melt point [2^] and triangle for PT pressure 0]) 
and other calculations (line plus small solid circles for MD 
[H and large circle for LBA+U Hp]). 



rite phase have been observed by calculating elastic con- 
stants. As a result, the melting temperature T m also in- 
creases linearly with pressure. Phonon dispersion results 
of the Fm3m phase are in good agreement with available 
experimental values. The LO-TO splitting at the T point 
is successfully reproduced by including the polarization 
effects. For the cotunnite phase, the imaginary mode 
along T—X direction and soft phonon mode along T—Z 
direction have been found at the equilibrium volume. 
The cotunnite to fluorite transition can be reached by 
firstly slipping the adjacent (100) planes relatively along 
the [001] direction to create a face-centered orthorhom- 
bic structure and secondly expanding the cell along the 
[010] direction and shrinking in the vertical directions 
to form the fee fluorite structure. Using QHA and the 
Debye model, we have calculated the Gibbs free energy, 
temperature dependences of lattice parameter and bulk 
modulus, entropy, specific heat, and P — T phase diagram 
of the important nuclear fuel material U02- 



IV. CONCLUSION 



In summary, we have carried out a first-principles 
DFT+ U exploration of the ground state properties as 
well as the high temperature/pressure behaviors of UO2. 
By choosing the Hubbard U parameter around 4 eV 
within the LDA+ U approach, the equilibrium state fea- 
tures for both the ambient Fm3m and the high-pressure 
Pnma phases of UO2 are shown to agree well with ex- 
periments. The Fmim^Pnma transition is predicted 
to occur at 40 GPa upon compression at K. At ambi- 
ent pressure, a transition temperature of 2069 K between 
the two solid structures is first obtained by the Debye 
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